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ABSTRACT 



A finite element shallow-water model is tested with two types of 
surface topography. The model uses rectangular subdivisions in a vorticity- 
divergence formulation, and a semi-implicit time discretization. In the first 
experiment an east-west ridge or valley is placed in a channel with east- 
west periodic conditions. Linear quasi-geostrophic solutions are derived 
with the rigid lid assumption. The Rossby waves are successfully simulated 
in the model with linear solutions as the initial conditions. The model phase 
speeds are very close to the analytic values when the latter are properly 
corrected. In the second experiment a ridge is placed across the channel and 
the Coriolis parameter is set to zero. The initial conditions consist of a 
uniform flow through the channel and constant free-surface height. The 
numerical simulations agree with hydraulic jump theory. In the jump cases 
the model predicts increasing wind speeds and decreasing free surface 
heights. Higher spatial resolution would be required to properly simulate 
the details of the hydraulic jump formation. 
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A. GENERAL 

Wilhelm Bjerkens (1904) was the first to point out that the future 
meteorological conditions can in principle be obtained by an integration of 
differential equations which govern the behavior of the atmosphere. Such 
an integration performed using numerical methods is called niLULS. Lic.Q.1 
w_ea t her jLLS. diQ.tL Q.ILi 

Richardson was the first to attempt a numerical weather prediction. 
After very long and time-consuming computations, he obtained a totally 
unacceptable result (Richardson, 1922). That wrong result, and 
Richardson’s estimate that 64,000 men are required to advance the 
calculations as fast as the weather itself is advancing, left some doubt 
about the practical use of the method. However, a number of developments 
that followed improved the situation. Mainly due to the work of Rossby in 
the late 1930' s, it became clear that even a rather simple equation, one that 
describes the conservation of absolute vorticity following the motion of air 
particles, suffices for an approximate description of large scale motions of 
the atmosphere. 

Finally, in 1945, the first electronic computer, ENIAC, was 
constructed. The absolute vorticity conservation equation, and ENIAC, 
were used by Charney, Fjortoft and Von Neumann in the late 1940' s for the 
first successful numerical forecast (Charney et al. 1950). Much faster 
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computers, and improved understanding of computational problems, now 
also enable long-term integrations of the basic primitive equations. With 
the introduction of each new generation of computers, the gap between 
numerical forecasts and atmospheric observations has decreased, but the 
rate at which this gap is decreasing appears to be leveling off. This 
indicates that technological improvements in computing power may not be 
the primary limitation to better numerical forecasts. 

During the past 15 years, there has been a significant effort within the 
numerical weather prediction area in developing limited- area, fine-mesh 
primitive equations models and applying them to operational, short range 
weather forecasts. An important practical motivation for the development of 
regional models has been the limited success of operational global models 
in the prediction of precipitation and severe weather. A parallel motivation 
for the development of regional models is their potential scientific value to 
researchers studying the structure and dynamics of mesoscale phenomena. 
(Keyser and Uccellini, 1987) 

B. BACKGROUND 

There are two fundamental methods of simulating the atmospheric 
flow; physical-models and mathematical-models techniques. With the 
physical-models technique, we construct scale model replicas of observed 
ground surface characteristics and insert them into a chamber such as a 
wind tunnel. The flow of air in this chamber is adjusted so as to best 
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represent the large scale, observed atmospheric conditions. Mathematical 
modeling, on the other hand, makes use of such basic analysis techniques 
as algebra and calculus to solve directly the equations governing 
atmospheric flows. 

Numerical solution of the equations of motion is performed using the 
grid point method for most applications. Following this method a set of 
points is introduced in the region of interest and dependent variables are 
initially defined and subsequently computed at these points. This set of 
points is called a grid. It is very important to note that most of the time the 
grid points are at fixed locations in the horizontal. This means that, 
according to the Eulerian system of equations, space and time coordinates 
are chosen as independent variables. 

With the grid point method, the most common way of solving the 
governing equations is to find approximate expressions for derivatives 
appearing in the equations. The required approximate expressions are 
defined using values of the dependent variables only defined at the grid 
points, and at discrete time intervals. Thus, they are formed using 
differences of dependent variables over finite space and time intervals ; that 
is the reason why this approach is called the finite difference method. The 
approximations for derivatives are then used to construct a system of 
algebraic equations that approximates the governing partial differential 
equations. This set of algebraic equations is to be solved, usually using an 
electronic computer, by a proper step-wise procedure in time. 
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A major limiting factor of finite difference approximations is the 
truncation error. That is, the smaller the grid interval, the smaller the 
truncation error. For a finite difference model to improve its accuracy, it 
would require increasing the grid matrix density. This would require 
increased computer storage and computational time. The problem here is 
not a simple one. It goes beyond numerical techniques and computer 
technology. For instance if we further reduce the grid spacing on the 7LPE 
(National Weather Service 7 Layer Primitive Equation Model) we do not get 
any significant improvement in the accuracy of the solution (Woodward, 
1981). The required additional computer capability can not be utilized using 
finite difference methods. Therefore, new numerical integration techniques 
must be investigated. 

Two alternative techniques, the spectral method and the finite element 
method, have been the subject of intensive research. Both the spectral and 
finite element methods require more computational time per forecast time 
step than does the finite difference method. That is why both the spectral 
and finite element methods must utilize efficient numerical techniques to be 
considered a viable option for numerical weather prediction. For long range 
weather predictions, the spectral method appears to be a natural method, 
when applied over the globe or hemisphere. This is due to the existence of 
efficient transforms for the nonlinear terms in spherical geometry. 
However, the spherical harmonics are globally rather than locally defined. 
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That is the reason why the finite element method appears to be more 
suitable for problems of more detailed limited area forecasting. 

The Galerkin Finite Element Method (GFEM) has the potential to 
increase efficiently the spatial resolution for the purpose of simulating 
accurately the small-scale processes. More specifically, the GFEM model 
used by Hinsman, 1983, has demonstrated desirable characteristics. It 
propagates atmospheric waves better than an equivalent finite difference 
model. It also allows variable-resolution grids and responds better than an 
equivalent finite difference model near the smallest gridlength. Moving 
grids can be achieved with no apparent noise generation. Finally, it can 
utilize direct solvers and is a natural choice for vectorization on large 
computers. Furthermore, the superior small-scale response of the GFEM 
indicates potential increase in skill for regional forecasting, and it truly is a 
viable option for simulation of atmospheric flow. 

C. OBJECTIVES 

The main objective of this study is to test the Galerkin Finite Element 
model which was developed by Hinsman (1983), with topography. This 
shallow-water model uses bilinear basis functions. An earlier study with 
the same model was carried out by Neta et al. (1986), in which the bottom 
topography changes linearly to the north. 

In this thesis two types of bottom topography will be considered in a 
channel domain encompassed by solid north-south walls with east-west 
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boundary conditions. In the first case the bottom is composed of two 
regions. These regions have a constant but opposite northward slope so 
that they can form either an east-west oriented ridge or an east-west 
oriented valley in the surface topography. Lines of constant bottom height 
run parallel to the x-axis, and pure geostrophic motion is possible only if 
the v-component of the velocity is identically zero. If we consider no mean 
flow and no beta effect, topographic Rossby waves can exist moving either 
east or west depending on the slope of the bottom. 

In the second case, a mean zonal flow passes over a ridge which 
extends north-south across the channel. It is clear now that the lines of 
constant bottom height run parallel to the y-axis. The Coriolis parameter 
for this case is set to zero, and the formation of hydraulic jumps is to be 
investigated for different values of the mean flow and peak of the 
topographic ridge. 
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A. HISTORICAL INTRODUCTION 

The finite element concept is to a large extent physical rather than 
abstract in nature, and it has been used in a variety of forms for centuries. 
The basic idea has always been to replace an actual problem by a simpler 
one. In other words the finite element method is the replacement of 
continuous functions by piecewise approximations, usually polynomials. 
Indeed the early geometers used "finite elements" to determine an 
approximate value of K. They did this by bounding a quadrant of a circle 
with inscribed and circumscribed polygons, the straight-line segments 
being the finite element approximations to an arc of the circle. In this way 
they were able to obtain extremely accurate estimates. Upper and lower 
bounds were obtained, and by taking an increasing number of elements, 
monotonic convergence to the exact solution would be expected. 
Archimedes used these ideas to determine areas of plane figures and 
volumes of solids, although of course he did not have a precise concept of 
a limiting procedure. The interesting point here is that while many 
problems of applied mathematics are posed in terms of differential 
equations, the finite element solution of such equations utilizes ideas which 
are in fact much older than those used to set up the equations initially. 

The modern use of finite elements really started in the field of 
structural engineering. Probably the first attempts were by Hrennikoff 
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(1941) and McHenry (1943) who developed analogies between actual 
discrete elements and the corresponding portions of a continuous solid. The 
term "finite element" was introduced later by Clough (1960) in a paper 
describing an application in plane elasticity. 

The engineers had put the finite element method on the map as a 
practical technique for solving their elasticity problems, and although a 
rigorous mathematical basis had not been developed, the next few years 
saw an expansion of the method to solve a large variety of structural 
problems. The workers in the early- 1960s soon turned their attention 
towards the solution of non-linear problems. Turner et al. (1960) showed 
how to use an incremental technique to solve geometrically non-linear 
problems, i. e., problems in which the strains remains small but the 
displacements are large. Stability analysis also comes into this category and 
was discussed by Martin (1965). Plasticity problems, involving non-linear 
material behavior were modelled at this time (Gallagher et al. 1962), and 
the method was also applied to the solution of problems in visco-elasticity 
(Zienkiewicz et al. 1968). 

Finally, besides the static analysis, dynamic problems were also being 
tackled, and Archer (1963) introduced the concept of the consistent mass 
matrix. Both vibration problems (Zienkiewicz et al. 1966) and transient 
problems (Koenig and Davids, 1969) were considered. Thus the period 
from its conception in early-1950s to the mid-1960s saw the method being 
applied extensively by the engineering community. Once it was realized 
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that the method could be interpreted in terms of variational methods, the 
mathematicians and engineers were brought together, and many extensions 
of the method to new areas soon followed. In particular it was realized that 
the concept of piecewise polynomial approximation offered a simple and 
efficient procedure for the application of the classical Rayleigh-Ritz 
method. The method had now become an important technique from both a 
practical and theoretical point of view, and the number of published paper 
using this method began to increase at a tremendous rate. 

In the area of meteorology also, the finite-element method has been 
successfully employed in the horizontal representation of atmospheric 
variables in numerical weather prediction and atmospheric modeling 
(Cullen, 1974a, 1974b, 1979; Hinsman, 1975, 1983; Staniforth and 
Mitchell, 1977, 1978). The finite element method when applied to 

meteorological equations gives very accurate phase propagation and also 
handles nonlinearities very well. 

B. FINITE ELEMENT APPROXIMATIONS 

In contrast to the finite difference schemes wherein the domain of 
interest is replaced by a set of discrete points, in the finite element method 
the domain is divided into subdomains called finite elements. The unknown 
function, let us name it u, is represented within each element by an 
interpolating polynomial which is continuous along with its derivatives to a 
specified order within the element. Generally, the interpolating function is 
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of lower-order continuity between elements than within an element. Thus 
the fundamental building block in the finite element method is the 
subdomain or element. 



C. THE METHOD OF WEIGHTED RESIDUALS 

There are several ways that can lead us to the same finite element 
formulation. A conceptually simple approach can be formulated using the 
method of weighted residuals. The two primarily special cases of the 
method of weighted residuals (MWR) are the Galerkin and the collocation 
methods. 

In the method of weighted residuals, the desired function u (•) is 
replaced by a finite series approximation as 



In general, the set of functions <})j (•), j = 1,2,...,N, can be defined over 
both the time and space domain and Uj, j = 1,2...,N, are undetermined 
coefficients. The equation (2-1) can be written 



N 



u (•) - a (•) = X u j 4>j <•)• 



( 2 - 1 ) 




N - 1 



( 2 - 2 ) 



j = i 
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where <j>j (•) satisfy the homogeneous boundary conditions. The functions 
<}>j (•) are chosen to be polynomials that satisfy certain of the boundary 
conditions imposed on the problem. These functions are variously denoted 
shape functions, basis functions, and interpolation functions, depending 
upon the discipline in which the method is being applied. Even if we 
choose the basis functions to satisfy all boundary conditions, they will not 
normally satisfy the PDE as well. If we now substitute the u (•) into the 
PDE, say Lu - f = 0, we can easily get 



where R (•) is a residual. 

Our objective at this point should be to select the undetermined 
coefficients Uj such that this residual is minimized in some sense. A 
straightforward scheme would be to set the integral of R (•) to zero as 



This scheme, however, generates only one equation for the N unknown 
coefficients Uj. It can be suitably modified by introducing weighting 
functions wj (•), i = 1,2, ...,N. 

Setting the integral of each weighted residual to zero yields n 
independent equations 



L u (•) - f = R (•), 



(2-3) 




(2-4) 



t V 
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1 1 R (•) w. (•) dv dt = 0, i = 1,2,...,N. 



(2-5) 



t V 



At this point we can solve equation (2-5), in theory at least, for N 
unknown coefficients. Equation (2-5) represents the general equation 
describing the MWR, and a multiplicity of schemes arise out of this one 
expression through the definition of the weighting functions w;. 

Among the MWR family of methods, the Galerkin, subdomain, and 
collocation schemes are most commonly encountered in practice. The one- 
dimensional weighting function for the Galerkin scheme is illustrated in 
Fig. 2.1, for the subdomain scheme in Fig. 2.2, and for the collocation 
scheme in Fig. 2.3. 

D. GALERKIN METHOD 

The Galerkin method results when the weighting function is chosen to 
be the basis function, as defined in (2-1). Thus we have 



The basis functions are formally required to be members of a complete set 
of functions. Because a complete set of functions can exactly represent 




( 2 - 6 ) 



t V 
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Fig. 2.1 Schematic representation of the one-dimensional 
weighting functions for the Galerkin method. (It is 
assumed that the chapeau function is used as a basis 
function). 
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Fig. 2.2 As in Fig. 2. 1 , except for the subdomain method. 
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Fig. 2.3 As in Fig. 2. 1 , except for the collocation method. 
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any function of a given class, the series of (2-1) is inherently capable of 
representing the exact solution as the number of terms in the series is 
increased. 

The requirement of completeness allows an alternative interpretation of 
the Galerkin formulation. A continuous function must be zero if it is 
orthogonal to every member of a complete set. The Galerkin method can be 
viewed as a scheme in which the residual is forced to zero in the sense that 
it is made orthogonal to the complete set of functions <{>i (•). 

In Fig. 2.1 the weighting function (and therefore the basis function) is 
a hatshaped, piecewise linear function, and because of its hatlike 
appearance, it is sometimes called a "chapeau" function. It is often 
encountered in the formulation of the finite element method. The chapeau 
function is the simplest of the basis functions in common use. 

E. TWO-DIMENSIONAL BASIS FUNCTIONS 

The extension of the weighted residual method to higher dimensions is 
relatively straightforward, provided that regular rectangular subspaces are 
employed. We can easily visualize the two-dimensional case using an 
elementary example, (Lapidus and Pinder, 1982). The discretized domain is 
simply illustrated in Fig. 2.4, and the two-dimensional basis function that 
is linear along each side in Fig. 2.5. 
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Fig. 2.4 Discretized domain for two-dimensional heat flow. 

Finite element nodes are indicated by small circles, and 
the elements themselves by Roman numerals. 
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Fig. 2.5 Two-dimensional basis function that is linear along 
each side. 
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Consider now the problem of two-dimensional, time-independent heat 
flow in a rectangular plate with a heat source located at the center of the 
plate. The governing equations for this problem will be given by 



£ ( T ) = T + T = Q 

v ' xx yy 


(2-7) 


T (x,2) = 1 


(2-8) 


T (0,y) = 1 


(2-9) 


T y (x,0) = 0 


(2-10) 


T x (2,y) = 0 


(2-11) 


Q (x,y) = Q w 5(x -1) 8(y - 1), 


(2-12) 



where x and y are Cartesian coordinates, Q w is the heat source, and 8 is the 
Dirac delta function. 

-L. 

Let us first define the trial functions: 

9 

T (x,y) = f (x,y) = X T j ( x »y)* (2-13) 

j= i 
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For convenience, it is advantageous to define basis functions in 
local (£,r|) coordinates where upon (2-13) becomes, for 

4 

T (x,y) = f (x,y) = X T j (2-14) 

j = i 

and <j)j (^ , T) ) are bilinear chapeau functions such as those illustrated in Fig. 
2.5. 



At this point we can easily formulate the integral equations using 
Galerkin's procedure. We can visualize the whole procedure as the 
requirement of orthogonality between the residual R and each basis 
function, as 

J R (x,y) 0. (x,y) dx dy = 0, i = 1,...,9, (2-15) 

A 



where 



R (x,y) = £ ( T ) - Q. (2-16) 

A is the domain of integration and £ is defined by (2-7). Substituting (2-7), 
and (2-16) into (2-15) yields 
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i = 



(2-17) 



J ( T xx + T yy - Q ) 4>i (X,y) dx dy = 0, 

A 

Applying Green's theorem to (2-17) to incorporate second- and third-type 
boundary conditions directly into the set of integral equations, equation (2- 

17) becomes 

J< T *^ + Vyi + Q ^ d y + + 

A S 

i = 1,...,9, (2-18) 

where l x and l y are direction cosines with respect to the normal to the curve 
S, the boundary of the domain A. In this problem, the second term of (2- 

18) will be used to conveniently define the zero-gradient Neumann 
boundary conditions of (2-10) and (2-11). 

If we now substitute the trial functions, as defined by (2-13) into 
(2-18) we can obtain the following set of algebraic equations 

J ( X T j (< t>xj ^>xi + ^yj ^yi ) + Q^i ) dx dy 

A J=1 

-J(f x l x +f y l y )(}) i ds=0, i = 1,...,9. (2-19) 

s 
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Because <j>i is defined such that it is nonzero only over elements adjacent to 
mode i (see Fig. 2.4 and 2.5) the integrations of (2-19) may be performed 
piecewise over each element and subsequently summed. Thus we can write 



4 9 

X J(X T j^xj ( J ) xi + ^yj^yi) + Q ( I > i) dxd y 

e=l Ae j-1 

-J (f x l x + f y l y ) <>. ds = 0, i = 1,...,9, (2-20) 

where A e is the area of element e, and S e is the curve bounding A e . 

It is now time to formulate the matrix equation. In general, the 
best way is to express (2-20) in terms of the local (^,tj) coordinate system 
to facilitate integration (see Fig. 2.6a, Fig. 2.6b). This is easily 
accomplished provided that the relationship between derivatives of <J>j in 
each coordinate system is readily available. To find this relationship, we 
have to employ the chain rule to obtain 
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Fig. 2.6a 



Rectangular element in (x,y) coordinates. 



RFPROODCED AT GOVERNMENT EXPENSE 




Fig. 2.6b Same element as Fig. 2.6a transformed into the 
(^ , tj ) local coordinate system. 
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3 X 3<t>j 3y 



3<|)j 3x 3<t>j dy 

drj dx 3^ 3 y 3^ 

This set of equations can be written as 



1 

i 




rn 


s 

iTTt 

1 




l 

< 

'w'' 

X 

1 


,F — TfW, 




1 






1 



or 



»j>, 

»jV 






( 2 - 21 ) 



where [J] is the Jacobian matrix. For our problem we can easily evaluate 
the [J] as 

0.5 0 . 0 ' 

0.0 0.5 

and the [J ]- 1 as 
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! 2.0 0.0 



0.0 2.0 



( 2 - 22 ) 



Equations (2-20) and (2-22) may be combined, which can give us 

4 + 1 + 1 Q 




- f (f 1 +f 1 )<|>.ds = 0, i = 

J v xx yy T i 



( 2 - 23 ) 



In changing the limits of integration we introduce the following 
relationship 



four integrals appearing in each element in (2-23) are identical for each 
element. The assembly process, that is, the transformation from element to 
global integrations, can be easily visualized using Table I which provides 
us with the relationship between Local Node Numbers and Global Node 
Numbers. 

The assembly procedure now calls for extracting information 
concerning the integration at the element level from the Table I. From 



dx dy = det [J] di; drj. 



(2-24) 



Because in our problem the elements are all of the same size, the 
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Element Nodes ' - — 




Global Nodes 
(Global Node Numbers) 




(Local Node Numbers) 


I 


II 


III 


IV 


1 


1 


2 


4 


4 


2 


4 


5 


7 


8 


3 


5 


6 


8 


9 


4 


2 


3 


5 


6 



Table I Relationship between Local Node Numbers and Global 
Node Numbers. 
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column I we can see that global nodes 4 and 5 correspond to element nodes 
2 and 3 in element I. In order to find the contribution to the global integral 
concerning these nodes, we have to create Table II which evaluates the 
integrals for one element in local coordinates using equation (2-23). From 
Table II the seeking contribution is found in row 2, column 3, and it is 
equal to -1/6. This value should be placed in matrix location (4,5) of the 
global matrix. However, this value is not final, because there is also 
information to be retrieved from element III. This information is located in 
the position (1,4) of the element matrix of Table II, and it is equals to -1/6. 
This value has to be summed with the previous integral value and the new 
value should to be placed in the same global matrix position (4,5). 
Combination of the two integral values yields the final value of -1/3, which 
can be seen in Table III in row 4, column 5. 

In general, the element coefficient matrix (Table III) is different 
for each element because of changes in either element geometry or 
parameter values. Moreover, it is often necessary to perform the 
integrations of (2-23) numerically. 
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1 


2 


3 


4 


1 


2/3 


- 1/6 


- 1/3 


- 1/6 


2 


- 1/6 


2/3 


- 1/6 


- 1/3 


3 


- 1/3 


- 1/6 


2/3 


- 1/6 


4 


- 1/6 


- 1/3 


- 1/6 


2/3 



Table II Evaluated Integrals for One Element in Local 
Coordinates for equation (2-20). 
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1 


2 


3 


4 


5 


6 


7 


8 


9 


1 

1 


2/3 


- 1/6 


0.0 


- 1/6 


- 1/3 


0.0 


0.0 


0.0 


0.0 


2 


- 1/6 


4/3 


- 1/6 


- 1/3 


- 1/3 


- 1/3 


0.0 


0.0 


0.0 


3 


0.0 


- 1/6 


2/3 


0.0 


- 1/3 


- 1/6 


0.0 


0.0 


0.0 


4 


- 1/6 


- 1/3 


0.0 


4/3 


- 1/3 


0.0 


- 1/6 


- 1/3 


0.0 


5 


- 1/3 


- 1/3 


- 1/3 


- 1/3 


8/3 


- 1/3 


- 1/3 


- 1/3 


- 1/3 


6 


0.0 


- 1/3 


- 1/6 


0.0 


- 1/3 


4/3 


0.0 


- 1/3 


- 1/6 


7 


0.0 


0.0 


0.0 


- 1/6 


- 1/3 


0.0 


2/3 


- 1/6 


0.0 


8 


0.0 


0.0 


0.0 


- 1/3 


- 1/3 


- 1/3 


- 1/6 


4/3 


- 1/6 


9 


0.0 


0.0 


0.0 


0.0 


- 1/3 


- 1/6 


0.0 


- 1/6 


2/3 



Table III Evaluated Integrals for the Four Elements of Fig. 2.4 
in Global Coordinates for Equation (2-20). 
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ii. Si£J2_IV 



We now solve the matrix equation (2-23) using Table III as 



2/3 


- 1/6 


0.0 


- 1/6 


- 1/3 0.0 


0.0 0.0 


0.0 




Ti 


- 1/6 


4/3 


- 1/6 


- 1/3 


- 1/3 - 1/3 


0.0 0.0 


0.0 




t 2 


0.0 


- 1/6 


2/3 


0.0 


- 1/3 - 1/6 


0.0 0.0 


0.0 




t 3 


- 1/6 


- 1/3 


0.0 


4/3 


- 1/3 0.0 


- 1/6 - 1/3 


0.0 




t 4 


- 1/3 


- 1/3 


- 1/3 


- 1/3 


8/3 - 1/3 


- 1/3 - 1/3 


- 1/3 


* 


t 5 


0.0 


- 1/3 


- 1/6 


0.0 


- 1/3 4/3 


0.0 - 1/3 


- 1/6 




t 6 


0.0 


0.0 


0.0 


- 1/6 


- 1/3 0.0 


2/3 - 1/6 


0.0 




t 7 


0.0 


0.0 


0.0 


- 1/3 


- 1/3 - 1/3 


- 1/6 4/3 


- 1/6 




Tg 


0.0 


0.0 


0.0 


0.0 


- 1/3 - 1/6 


0.0 - 1/6 


2/3 




t 9 
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f (f 1 +f 1 ) <t>, ds 
j v xx yy T i 

s 

s 

f ( t 1 + 1 1 ) <J>, ds 

J v xx y y 7 T 3 



■Qw 



f ( t 1 + f 1 ) <t>. ds 

J y xx yy /Y 6 



o.o 



( 2 - 25 ) 



f ( t 1 + f 1 ) <> 0 ds 

J ' xx yy /r 9 
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The line integrals appearing on the right-hand side of (2-25) 
represent flux boundary conditions. In general, when Dirichlet boundaries 
are specified, it is necessary to expand and evaluate the line integrals. 
However, when we use functions, that have continuous derivatives up to 
the Oth order as bases, one can condense rows containing the known 
temperature values out of the matrix equation. The reduced matrix equation 
is 



W ‘ (b) = (f). 



where 




4/3 


-1/3 


-1/6 


-1/3 


-1/3 


8/3 


-1/3 


-1/3 


-1/6 


-1/3 


2/3 


-1/6 


-1/3 


-1/3 


-1/6 


4/3 



{b} = 




(6-26) 



(6-27) 



( 6 - 28 ) 
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{f} = 



1/2 

-Q +5/3 

0.0 

1/2 



( 6 - 29 ) 



Solving (2-26) directly we can easily get 



— 

t 4 




0.783 


t 5 




0.525 


t 7 




0.655 


t 8 




0.783 



( 2 - 30 ) 



The coefficients T 4 , T 5 , T 7 , and T 8 represent the values of the 
temperature at nodes 4, 5, 7, and 8 because the basis functions are defined 
such that 4 >i is unity at node i and zero elsewhere. 

At the same time the analytical solution for our problem is 



T = U + V + W, 



(2-31) 
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where 



2 y sin[(n + l/2)7t(a - y)/a] cosh[(n + l/2)rc(a - x)/a] 
^ n = 0 (n + 1/2) cosh(n + 1/2)7C 



2 Y cosh[(n + l/2)7ty / a] sin[(n +l/2)7tx / a] 
K n = o (n + 1/2) cosh(n + l/2)Jt 



( 2 - 33 ) 



w _2 y sin[(n+l/2)7rx/a]sin[(n+l/2)TC^/a]sinh[(n+l/2)7c(a-y)]cosh[(n+l/2)7rri/a] 

^ n = o (n + 1/2) cosh[(n + 1/2 )tc] 

(2-34) 

where a is the length of the side of the square (in our case equals to 2), 
and tj are the x and y locations of the heat source, respectively. The 
analytical solution yields 



a 




' 0 . 756 ' 


T 2 




0.127 


T S 




0.719 


a 




0.756 



(2-35) 
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It is obvious that the numerical solution at the singular point (1,1), which 
corresponds to the nodal location 5, is not very accurate. This is not odd 
since we are attempting to represent a rapidly varying function with only 
four bilinear elements. As we mentioned before the series of (2-1) is 
inherently capable of representing the exact solution as the number of terms 
in the series is increased. 
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Ill, THE SHALLQ^^ATER_MQDEL-AriD-TH£_TQPQIIRAPHI£ 

ROSSBY WAVE 



A. GENERAL 

In order to study motions of atmospheric and oceanic relevance, we 
can use a shallow, rotating layer of homogeneous fluid which is 
incompressible, and inviscid. This model (Shallow-Water Model) ignores 
completely the presence of stratification, but experience has shown that it 
is capable of describing important aspects of atmospheric and oceanic 
motions. Because of this, it is useful to deal with shallow fluid systems in 
trying to determine general principles of hydrodynamical behavior of the 
atmosphere. 

The major physical characteristics concerning the shallow-water 
model, are found to apply well for more complex systems. If we were to 
analyze the types of wave motions associated with more complex forms of 
the primitive atmospheric equations, (e.g., with vertical stratification, 
compressibility, etc.), we would find essentially the same types, namely 
gravity-inertia waves and Rossby waves, and deep and shallow motions 
would exist simultaneously. 

B. THE SHALLOW- WATER MODEL 

Let us consider a sheet of fluid with constant and uniform density as 
illustrated in Fig. 3.1. The height of the surface of the fluid is given by 
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Fig. 3.1 The Shallow-Water Model. 
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h(x,y,t). We also model the body force arising from the potential 4> as a 
vector, g, normal to the z = 0 surface. The z-axis coincides with the 
rotation axis of the fluid, so that in this particular case the Coriolis 
parameter f is simply given by 2Q. The rigid bottom topography is defined 
by the surface z = he(x,y), which is not a function of time (t). Finally, we 
assume that the fluid is inviscid ([1=0), that is, only motions for which 
viscosity is unimportant are considered. 

Also we suppose that a characteristic value for the depth can be 
sensibly chosen, say D, and D also characterizes the vertical scale of the 
motion as well. In the same way we consider that a characteristic horizontal 
length scale for the motion exists, which we call L. The fundamental 
relationship which characterizes the shallow-water theory is given by 

5 = ^ « 1. (3-1) 

Here, it is important to note that the fluid is rotating, so that Coriolis 
accelerations can be important. The fluid layer is flat rather than forming a 
spherical shell, and its major physical deficiency as mentioned above is the 
absence of the density stratification which is present in the real 
atmosphere. 
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C. SMALL. AMPLITUDE MOTIONS 



Dealing with the Shallow-Water Model, and since by hypothesis 5<<1, 
we are able to express the hydrostatic approximation in the following form: 



3p 

dz 



= -pg- 



(3-2) 



We can integrate equation (3.2) resulting in 

p = - p g z + A(x, y, t). ( 3 - 3 ) 

Applying now the obvious boundary condition 



p(x,y,h) = p Q , 



(3-4) 



where po is a constant, we can easily get 
p = pg(h-z) + p fl . 



(3-5) 



At this stage, we can clearly see that the horizontal pressure gradient 
is independent of z, i.e. 



3p 3h 

aT =pg 



(3-6) 
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Op Oh 

o7 _pg o7’ 



(3-7) 



in such way that the horizontal accelerations must be independent of z. That 
is why we can also assume that the horizontal velocities themselves remain 
independent of z, if they are so initially. Applying now the Taylor- 
Proudman theorem we are able to write the horizontal momentum equation 
as 



Ou 

dF 



+ 



Ou Ou „ 

U -=r— +V-=— - fv = - 
dx dy 



dh 

g dx ’ 



(3-8) 



dv dv dv „ dh 

dT +u d7 +v d7 +fu = ' g d7' 



(3-9) 



The specification of incompressibility for the Shallow-Water Model, 
decouples the dynamics from the thermodynamics and reduces the equation 
of mass conservation to the condition of incompressibility given in the 
following form: 



du dv dw 
dx dy dz 



(3-10) 



Having in mind that u and v are independent of z, equation (3-10) can be 
integrated in z as 



41 




(3-11) 



Equation (3-11) can be further manipulated using the condition of no 
normal flow at the rigid surface z = hg resulting in 



Our next step now should be to linearize the set of equations (3-8), (3- 
9), and (3-15), by studying small-amplitude motions. This is very 
important because the presence of solutions representing free oscillations 
often demonstrates fundamental mechanisms which occur in more 
complicated situations as we mentioned before. 

Let us now consider the thickness of the fluid layer in absence of 
motion to be given as 



lr + lr l(h ' h B )ul % ,(h ' h B )v ) = °- 



(3-12) 



If we now consider the total depth (H) is given by 



H = h - h, 



(3-13) 



then the equation of mass conservation (3-12) becomes 




(3-15) 
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H(x,y,t) = H 0 (x,y) + r|(x,y,t), 



(3-16) 



where rj << Hq. Furthermore 



3 u h 

3t >>U H- VU H- 


(3-17) 



or in other words u and v are considered small enough. Linearizing now the 
equations (3-8), (3-9), and (3-15), we can obtain 



9u 3r| 

3T- fv = - g 3T' 


(3-18) 


9v dr) 

3T + fu = - g 37- 


(3-19) 




( 3 - 20 ) 



where all the quadratic terms in the dynamical variables u, v, "q with 
respect to the linear terms are ignored. If we now define the linearized 
mass flux vector given by 



U = iU+jV, 


(3-21) 


where 
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( 3 - 22 ) 



u = u H 0 , 


( 3 - 22 ) 


V = vH 0 , 


( 3 - 23 ) 


equations (3-18), (3-19), and (3-20) become 

IT 


( 3 -24 ) 




( 3 -25 ) 


9 ri au av „ 

3 — + - 5 — - 5 — = 0 . 
at ox dy 


(3-26) 



At this point it is possible for us to obtain an equation in the single 
variable r\ as 



a a 2 

^ [(— + f 2 ) 0 - V . (C 2 Vti)] - g f J (H 0 ,ti) = 0 , 
at 


( 3 - 27 ) 


where 

C 2 „=gH 0 , 


( 3 - 28 ) 
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( 3 - 29 ) 




It can be seen that the Jacobian term is the effect of the geostrophic 
wind blowing across isobaths. For low frequency motion and small 
variation in Ho, the first term in (3-27) represents the time rate of change 
of the quasigeostrophic vorticity. The motion that results from this kind of 
balance is a topographic Rossby wave which we examine in more detail in 
the next sections. The velocities components u and v can be also found in 
terms of T|, given as 



D. THE TOPOGRAPHIC ROSSBY WAVE 

Following Pedlosky (1987), let us consider that Ho in equation (3-27) 
varies slightly in the y-direction given by 




( 3 - 30 ) 




(3-31) 




( 3 - 32 ) 



where D is a constant. 
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s is the slope of Ho in the y-direction, and 
L is the width of the channel. 

In this particular case (Fig. 3.2), pure geostrophic motion is possible 
only if v is zero. Also, lines of constant Ho are parallel to the x-axis. 
Motion across the isobaths of fluid columns will cause them to stretch or 
contract. Therefore, there is a possibility of a different mode of motion to 
exist depending on the combined effect of rotation and bottom slope. 

Assuming 



s « 1, 



( 3 - 33 ) 



and T| to be of the following form 



T| = Re { T|*(y) exp[i(kx - crt)] } , 



(3-34) 



substitution in (3-27) yields 




9 V f S 

- k (1 - s— ^ k] =0, 




( 3 - 35 ) 



with the following boundary conditions: 
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Fig. 3.2 The infinite channel of width L, rotating with angular 
velocity f/2. 
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( 3 - 36 ) 



dri fk . 

-~ + — ri =0, 

d y a 



on y = 0, L. 

If we now assume that 






( 3 - 37 ) 



which appears to be a very realistic approximation, then (3-35) becomes 



•*L +[ 2'^. k > fi k ,V = o. 

dy 2 L dy cl L o 



( 3 - 38 ) 



Here only in terms where s can be compared with quantities of order 
unity has it been neglected. Solving (3-38) we can get 



11 = ex P(-|f) t A sin(ay) + B cos(ay)]. 



( 3 - 39 ) 



where a is given by 



a 2 - f 2 



.,2 s v fks 

a = / (k +— ) 

S> 4L CtL 



(3-40) 
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Applying now the boundary conditions given by (3-36) we can get the 
following eigenvalue problem 



(a 2 - f 2 ) (a 2 - k 2 d 0 ) sin (aL) = 0. 



(3-41) 



It is obvious, that to the lowest order in s, the slope does not alter the 
Kelvin mode (that is, because the factors multiplying sin (aL) are the same 
as for the case of a flat bottom). The roots corresponding to the zeros of 
sin (aL) are given by 



fksd 



2 2 



f 2 



. d (k 2 + = o. 

La 0 i . 2 r 2 



( 3 - 42 ) 



There are two separate classes of solutions to (3-42) 

a. The first class has frequencies each of which exceeds f. In this 
class the term in s is negligible, and to O(s) we obtain the Poincare modes, 
i.e.. 



2 2 

o = f 2 + C Q (k 4 — ) 4- O(s), n = 1,2,.... ( 3 - 43 ) 

L 



Equation (3-43) shows that the high-frequency Poincare waves are also 
essentially unaffected by the small bottom slope. 
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b. The class having frequencies o = O(s), for which the first term 
in (3-42) is negligible, while the second is of 0(1). 

The other solution leads us to the dispersion relation for the 
topographic Rossby wave, i.e.. 



o = 



- 4 > 



2 2 
,2 n 7t 
k + + 



r 



n = 1,2,..., 



(3-44) 



obtained by neglecting the first term of (3-41). 

The maximum Rossby-wave frequency should be given by setting 



nV f\l/2 

k = k „ = (— + — > • 
L 2 C$ 



(3-45) 



for which we get 



o = a = 

max 



2 2 f 2 L 2 1/2 
[n 7t + •] 



'0 



(3-46) 



From equation (3-46) it is clearly seen that the Rossby-wave frequency 
is always less than f. Another important feature of the Rossby wave is that 
its phase speed in the x-direction, which is given by 
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s f 

T~ 



( 3 - 47 ) 



is negative for s > 0, and positive for s < 0. Therefore, the wave 
propagates (in the Northern Hemisphere) parallel to the topography, with 
the shallowest water on its right. Also, for high wave number, i.e., small 
scale, the frequency a decreases with increasing wave number as shown in 
Fig. 3.3. 

The dynamical fields for the Rossby-wave to lowest order are given by 



ntty 

tj = tJq sm(— =-£-) cos(kx - at + (J>) + O(s), 

L 



(3-48) 



u = " f^TT^ 0 cos (“^ cos(kx - at + (J>) + O(s), 



(3-49) 



g n7Ty 

v = - y sin(kx - at + <J>) + O(s), 



( 3 - 50 ) 



where small terms of O(s) have been neglected consistently. 
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Fig. 3.3 A schematic representation of the dispersion diagram 
for Poincare, Kelvin, and topographic Rossby waves 
in a channel. 
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From equations (3-48), (3-49), and (3-50), it follows that to lowest 
order in s, and therefore o/f, the fields of motion in the Rossby wave 



satisfy 




(3-51) 




( 3 - 52 ) 



which are the time-independent forms of (3-31) and (3-32). Because of the 
close relationship between rj, and the pressure field, equations (3-51) and 
(3-52) represent the geostrophic relation for the horizontal motions. 



plays the role of the Rossby number here, the velocity fields, though 
changing with time, remain continuously in geostrophic balance with the 
pressure field. However, the flow is not exactly geostrophic, for then the 
flow would be restricted to travel parallel to the isobaths, i.e., v would 
vanish. Even though the velocity fields are geostrophic to lowest order, it 
is the very small departures from geostrophy that give rise to the wave. It 
is the small cross-isobath flow, which is a nongeostrophic effect, which 
produces the oscillation. The Rossby-wave, whose existence requires both 
s and f to be non zero, is a low-frequency wave oscillation; its period is 
greater than a rotation period. 



At this stage, it is clearly seen that to lowest order in a/f, which 
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IV_HXDRAiLLI£.l]IM£S._IN_Ri2TAlINii 

AND NON ROTATING SYSTEM S 



A. GENERAL 

A strong and relatively warm wind known as chinook occurs from time 
to time in areas along the eastern slope of the Rocky Mountains. Similar 
phenomena are often observed in the Owens valley on the eastern side of 
the Sierra Nevada. Also cold fronts, approaching the area of Alps, moving 
from north or west, many times undergo severe deformation. This 
deformation may result in a number of important weather events including 
blocking and splitting of air flow on the upstream side or even triggering of 
lee cyclogenesis in the downstream flow. 

Tepper (1952) has proposed that squall lines are modified hydraulic 
jumps, but the idea that there is a certain link between downslope winds 
and hydraulic jumps was proposed by Long (1953). Long suggested that 
the mechanism that produces the hydraulic jump may be similar to one that 
produces strong waves and downslope winds in the atmosphere. Houghton 
and Kasahara (1968) have also proposed other mesoscale hydraulic 
analogies. Williams and Hori (1970) observed a delay in the formation of 
hydraulic jumps in case the Rossby number was less than 0. 1 in a rotating 
system. 

However, it has been difficult to confirm this hypothesis about jump 
formation, because there are significant differences between the atmosphere 
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and the simple fluid systems used in the hydraulic theory. One of the major 
reasons for this uncertainty appears to be the vertical propagation of the 
wave energy, which may occur in the real atmosphere but not in fluids 
bounded by a rigid or a free surface. 



B. HYDRAULIC JUMPS IN A NONROTATING SYSTEM 

Let us first consider the one-dimensional Shallow-Water Model. Let us 
also consider a mean flow over an isolated rigid orographic ridge as shown 
in Fig. 4. 1. The horizontal momentum equation is given by 



3u 3u 

^— + Utt— 

3t 3x 



+ SaT< h + h M> =0 - 



(4-1) 



The continuity equation also, is given by 




3h 3u 

37 + h 37 



= 0 , 



(4-2) 



where h denotes the depth of the fluid, and 

hM is the height of the rigid ridge, which is a function of x. 

At time t = 0, the fluid is set in motion from rest so that for infinite x, 
it has a constant zonal flow uo- After sufficient time has elapsed, the 
solution in the neighborhood of the rigid ridge would be given by the 
steady state solutions of equations (4-1) and (4-2). 
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HCPROmJCED AT GOVERNMENT EXPENSE 




X ► 



Fig. 



4. 1 A cross section view of the one-dimensional Shallow- 



Water Model with a mean flow u. 
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Following Houghton and Kasahara (1968), the steady state solutions for 
the variables u and h are given by 



u 2 + 2gh + 2gh M = C 1 , (4-3) 

hu = C 2 , (4-4) 

where Ci and C 2 are constants. If a flow without a hydraulic jump is 
considered, both Ci and C 2 are determined by the velocity uo, and the 
height ho of the flow in the region far from the ridge, so that 



C i = u o + 2gh 0 , 


(4-5) 


^2 = h o u o- 


(4-6) 



At this stage we define the following dimensionless parameters Fq, F, 



R, and U* as 




•n 

0 

11 

irtTl 

EX 0 
O 


(4-7) 


2 F o 
F -f+1. 


(4-8) 
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(4-9) 



R = 



‘M 




U, 



u 




(4-10) 



If we eliminate h from equation (4-3) by using equation (4-4), we are able 
to obtain 



(F 2 - 1)U 3 , + (R-F 2 ) U* + 1 = 0. 



(4-11) 



Assuming that the mean flow uo is always positive, then F > 1. Dividing 
equation (4-11) by (F 2 - 1) we can get 

R F 2 1 

U* + (-— ) U, + — =0. ( 4 - 12 ) 

K - 1 F - 1 

Equation (4-12) is a very important relation among the variables U*, 
R, and F. We can easily plot the solution U» to equation (4-12) for given 
values of F 2 . For instance if F 2 equals to 1.02 (Fo = 0.2), equation (4-12) 
becomes 

U 3 + (50R-51)U, + 50=0. (4-13) 
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We can consider equation (4-13) as determining the values of R that are 
possible with a given value of U*, as shown in Fig. 4.2. It is interesting 
here to note the singular behavior near R = 0. In case of F = 1.125 (Fo = 
0.5), equation (4-12) becomes 

U* + (8R - 9) U* + 8 = 0, ( 4 - 14 ) 

which leads to Fig. 4.3. Also, if F = 3.0 (Fo = 2.0), equation (4-12) can 
be written in the following simple form 

U* + (0.5R - 1.5) U* + 0.5 = 0, ( 4 - 15 ) 



which gives Fig. 4.4. 

It is very important here to note, that if R has values lower than 
Rcriticab where R cr iticai is given by 



R 



critical 




1.8899 (F 2 - 1) 1/3 , 



(4-16) 



then three real roots exist, as we clearly see in Figs. 4.2, 4.3, and 4.4. In 
the case where R is greater than R C riticab no physically meaningful solution 
exists (no real roots). In other words, in this particular case we expect the 
occurrence of a hydraulic jump. 
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REPRODUCED AT GOVERNMENT EXPENSE 




Fig. 4.2 Parameter R as a function of solution U*, as given by 
equation (4-12) for F 2 = 1.02. 
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Fig. 4. 3 As in Fig. 4.2, except for F 2 = 1. 125. 
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Fig. 4.4 As in Fig. 4.2, except for F 2 = 3.0. 
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If we plot equation (4-16), as shown in Fig. 4.5, we are able to find three 
distinct domains. In both domains I and III, the parameter R (function of 
x), has values lower than R C riticai- In this case, a real solution of equation 
(4-16) there exists, which is physically meaningful. On the other hand in 
domain II, the parameter R is greater than R cr iticai> so no physical solution 
exists. 



C. HYDRAULIC JUMPS IN A ROTATING SYSTEM 

Although none of the numerical experiments concerning hydraulic 
jumps (results section), involves rotation, it is useful to explore the 
particular effects of rotation on the formation of hydraulic jumps. The basic 
equations for a homogeneous, one-layer, inviscid fluid (Williams and Hori, 
1970), are given by 



3u 3u 3h 

3r +u & +s 3r fv=0 ’ 



(4-17) 



3v 3v 

TT— + U -=r— + fu = 0, 
3t 3x 



(4-18) 




3h 3u 

37 +h 37 



= 0 , 



(4-19) 



where h is the depth of the fluid. 
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Fig. 4.5 Classification of asymptotic flow conditions as a 
function of the maximum height of the topographic 
ridge, R critica ], and initial flow speed parameter F. 
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We perform a scale analysis by expressing the independent variables 
as 

t = T t\ ( 4 - 20 ) 

x = L x'. (4-21) 

The dependent variables are also broken up and scaled as 



u = U u', 


(4-22) 


v = V v\ 


(4-23) 


/TT 




h = h + / — h\ 

m V g 


(4-24) 


where h m represents the mean depth of the fluid. 


Also, the Rossby number 


and the Froude number are defined as 




R 0" fL 


(4-25) 


and 






(4-26) 
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In order to include the results obtained by the previous analysis which 
are valid for a nonrotating system, we examine the case where Fo ^ 1 and 



Fq ^ Rq* The appropriate time scale for doing this is 



T = L 

(s h J 1/2 


(4-27) 


and the appropriate v scale 
F o 

V-jpU. 

R o 


(4-28) 



The nondimensional equations for this case are given by 



0u’ _ ,9u’ 9h’ Fo , n 

+ F n u + TT— v' = 0, 

3t 0 3x dx 0 2 

R o 


(4-29) 


3v’ 3v* 

3f + F ° u ' ax' + 


(4-30) 


3 3 %F o( u'f + h'^: )+ ^:=o. 


(4-31) 



It is obvious that hydraulic jumps can be formed through the action of the 
nonlinear terms in equations (4-29) and (4-31). Even when Fo is small 
(<<1), they can produce a jump in a nonrotating system. If the Coriolis 
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term in equation (4-29) dominates then a hydraulic jump may be prevented. 
This implies that 



The form of the curve dividing the jump region from the nonjump 
region should be given by 



where A is a constant. The numerical solutions show that the range for A is 
from 6.0 to 7.5, which appears to be in agreement with Houghton's 
analytical curve corresponding to A = 6.5. For Fo ~ 1, the above scale 
analysis does not apply, since Fo is then greater than Ro. 

When both Fo ~ 1, and Ro ~ 1, all terms in the equations (4-29), (4- 
30), and (4-31) are of the same order. In this case, jumps are expected to 
form. On the other hand, if Ro << 1, the proper time scale is 1/f, while the 
proper v scale should be V = U. The nondimensional equations can be then 
rewritten as 




(4-32) 




(4-33) 




(4-34) 
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3v' 3v' 

aT + R ° u ’ d7 + u ’ = °’ 



( 4 - 35 ) 



9u' 



9u’ 



9u’ 1 9u’ 



, + RrJu' ■=; — r + h' -=r — 7 + .. 

9t 0 3x 3x F n 9x 



-]-o. 



( 4 - 36 ) 



If we now neglect all terms in Ro, the resulting equations describe an 
inertial oscillation in u and v. We do not expect a hydraulic jump to form in 
this case except perhaps after a very long time. 
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V. MOD EL D ESCR IPTI ON 



A. GENERAL 

There are two possible choices of increasing resolution, where desired 
or required, concerning a triangular subdivision. We can use near- 
equilateral triangles (Cullen, 1974b) or equilateral triangles (Hinsman, 
1975). Both have the advantage of almost perfect wave propagation 
characteristics. For the same problem we can also use rectangular 
subdivisions, although it is not obvious which subdivision is most 
suitable. However, a major advantage of the rectangular subdivision 
(shown in Fig. 5.1) is that it allows algorithms to be developed which take 
full advantage of vector processors. The interesting point here is that the 
Galerkin method has to utilize efficient numerical techniques to be 
considered as a viable option for numerical weather prediction. 

There are several solution procedures available for the Galerkin 
method. No particular attempt is made to optimize computational efficiency 
in this research model. A direct solver is employed using a Gaussian 
elimination procedure. The matrices from the Galerkin procedure are 
decomposed into upper and lower block tri-diagonal matrices. A 
preprocessing, representing the forward substitution stage, can be done 
once. Any time a solution is desired, a back substitution has to be 
performed. That is why the required coefficients for the backward step 
must be stored in a very efficient manner. This particular algorithm 
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REPRODUCED AT GOVERNMENT EXPENSE 



Fig. 5.1 Rectangular uniform subdivision 
Cartesian coordinates. 



for a channel in 
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represents a 'skyline' solver, referring to the compact method of storing 
only those coefficients required. This method has both the desired level of 
accuracy and a high degree of computational efficiency. 



B. EQUATION FORMULATION 

In order to integrate the equations governing the free-surface height 
and velocity of an inviscid hydrostatic incompressible fluid we can write 



3 (]> 

dt 



+ u 




3<J) 



, . du 3v 

, " ( ar + a7 )=0 



( 5 - 1 ) 



3u 

It 



+ u 






( 5 - 2 ) 



dv dv dv d<b 

TT— + U •=— + V “ + fU + -=— = 0 

dt ox d y oy 



( 5 - 3 ) 



where <() is the geopotential height, 

u is the east/west component of the wind, 
v is the north/south component of the wind, and 
f is the coriolis parameter. 

We now assume that the geopotential height <j>, in absence of motion is 
<3>. Then, in general 
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<J) ( x, y, t ) = 0 ( x, y ) + <j>' ( X, y, t ), 



(5-4) 



where <t> is the mean, and 




<J)' is a perturbation from the mean. 


The governing equations now 


can be written 




4>' t + d>D + (u<J>') x + (v<}>') y = 0, 


(5-5) 


u ( + <J)' x + K x - vQ = 0, 


(5-6) 


V i + ^'y + K y + uQ = 0’ 


(5-7) 


where 




„ 3u 3v 

d= 5T + 37’ 


(5-8) 


1 2 2 
K = - (u + v ), 


(5-9) 


dvdu 

Q ~dx~ dy + f ' 


(5-10) 



Because of the rapidly moving gravity waves, the stability condition 
for a numerical integration normally requires a much smaller time step than 
for the simple advection equation since d> 1/2 >> U. Similar results may be 
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expected with the the more complete equations actually used in numerical 
weather prediction. Since the gravity waves are usually relatively 
unimportant in large-scale weather forecasting, the small time step required 
for computational stability increases the computing time considerably with 
little or no compensation by way of increased accuracy, perhaps even a 
loss. On the other hand, implicit differencing schemes, which may have no 
restriction on the size of the time step, have the serious disadvantage of 
requiring the inversion of a large matrix. 

A semi-implicit scheme has the great advantage of permitting a 
relatively large time step without unduly increasing computation time. In 
other words the semi-implicit scheme slows artificially the propagation 
speed of the fastest gravity waves, which allows a much larger time step 
than the normal Courant-Fredrich-Lewy (CFL) stability criterion and also 
offsets some of the extra computational expense required to solve the 
system of equations assembled at each step. For this reason it is now 
necessary for the implementation of the semi-implicit time discretization to 
rewrite our equations in terms of a velocity potential % and a 
streamfunction \| / defined by 

u =X x -V y , (5-11) 

V = *y + V x , (5-12) 
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with the following diagnostic relations 




(5-13) 




(5-14) 



In other words, we can use vorticity and divergence or velocity 
potential and streamfunction as variables instead of velocity components. 
This means that second-order derivatives appear in the equations and 
Poisson equations have to be solved. Using linear elements, the scheme 
obtained for the Poisson equation is very similar to the finite difference 
scheme and can be inverted by the same technique. Cullen and Hall (1979) 
showed that the accuracy of the Galerkin finite element method solution 
was more accurate for the vorticity-divergence formulation of the shallow- 
water equations than for an increase in resolution with the primitive 
formulation. This unstaggered vorticity-divergence together with staggered 
variable formulation gives the best treatment of geostrophic adjustment for 
small-scale features (Williams and Schoenstadt, 1980). 

Following the vorticity-divergence approach and dropping the primes 
for the rest of the chapter, the equations become 



<t> t + <t>D + (u(p) x + (v(J))y = 0 , 



(5-15) 
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(5-16) 



C t + V 2 <|> + (uQ) x + (vQ) y = 0, 

D t + V 2 <t) + V 2 K-(vQ) x + (uQ) y = 0. (5-17) 



where V 2 is the Laplacian operator and C, is the relative vorticity given by 



3v du 



(5-18) 



Now we express the velocity as the sum of the rotational and 
irrotational components 



V = V + V 

V X ’ 



(5-19) 



where 



V = k x Vw, 

v 



( 5 - 20 ) 



V =Vy 

x 



(5-21) 



and then we are able to rewrite the equations using 



9u dv 
dx + dy 




(5-22) 
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r 3v 9u _ 2 


(5-23) 


obtaining the following 




<j) t + ov 2 X = - (u0) x - (v0) y 5 


(5-24) 


(vV) t = -(uQ) x -(vQ) y> 


(5-25) 


(V 2 % ) t + V 2 0 = (vQ) x - (uQ) y - (K) x - (K) y . 


(5-26) 



Manipulating now the last equation (5-26) and taking care of the 
bottom topography (assumed to be not a function of time), we can easily 
obtain 



0 t + <J>V \ = - [u (0 - 0 B )] x - [v (0 - 0 B )] y ^ 


(5-21) 


(V 2 y) t = -(uQ) x -(vQ) yi 


(5-28) 


V 2 (% l + 0) = (vQ-K x ) x -(uQ + K y ), 


( 5 - 29 ) 



where 0 b is the bottom topography defined by the rigid surface z = he (x,y) 
not a function of time (t). 
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We now define the domain of integration to be a channel encompassed 
by solid north-south walls with east-west cyclic boundary conditions. The 
boundary condition at the walls should be 



where n is the outward pointing normal vector. 

Here it is interesting to note that, in rewriting the equations in this 
form, we have increased their order in x and in y from first to second order 
and should expect that it may be necessary to impose further boundary 
conditions. However, since we have sufficient boundary conditions, 
already, any further specification must not be arbitrary but should be a 
consequence of the previous formulations. Along the walls, the v 
component of the velocity has to be equal to zero, resulting in 



The zonal and meridional components of the wind can be written now as 



V . n = 0, 



( 5 - 30 ) 



<)> = - fu . 
Y y 



(5-31) 




(5-32) 




( 5 - 33 ) 
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Since the v component of the velocity has to be zero along the north 



and south walls then the obvious boundary condition should be 



V x + X y = 0, 



(5-34) 



and we can satisfy the above condition by simply setting 

V = 0, (5-35) 

when solving the vorticity equation, and 

X y = 0. (5-36) 

when solving the divergence equation. As a matter of fact this is an 
overspecification but equation (5-34) would be difficult to apply. Our 
initial conditions of course, must be specifically selected to satisfy both 
equations (5-35) and (5-36). 

As we mentioned before we use a semi-implicit time discretization 
scheme for reasons of computational efficiency. Basically, this semi- 
implicit scheme is simply a modified leapfrog scheme giving a net saving in 
the computational time required to make a forecast for a given time. The 
way in which this is accomplished is to evaluate certain terms implicitly as 
a mean over times (t - At) and (t + At) rather than at time t. 
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Following this approach we evaluate all the terms on the left hand side 
of the equations as an average at times (t - At) and (t + At), and all the right 
hand side at time t. The prognostic equations then become 

4> t + <D[V \ (t - At) + V \ (t + At)] = - [u (<]> - <|> B )] x - [v (<]) - ()) B )] y , ( 5 - 37 ) 

V 2 (\|/ t ) = -(uQ) x -(vQ) y) (5-38) 

V 2 [% t + <J> (t + At) + <J) (t - At)] = [(vQ) - K x ] x - [(uQ) + K y ] y . ( 5 - 39 ) 



If we now solve equation (5-39) for % (t + At) and substitute into equation 
(5-37) we can finally get the following set of equations 



V V r = [(vQ) - K ] - [(uQ) + KJ 



O (At) 



X J X 



(J> (t - At) 



y J y . .. .2 
<E> (At) 



V y(t - At) 1 

+ — - + {[U (9 - <I> B )] X + [V (0 - Vly} , 

At 0 At y 



(5-40) 



2 , * 



v «|> +X t ) = [(vQ)-K x ] x -[(uQ) + K y ] y> 



(5-41) 



V (v t ) = - (uQ) x - (vQ) 



(5-42) 



where <))* is given by 
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$ = [ 0 (t - At)-0(t + At)]/2. 



(5-43) 



At this stage our initial system of three equations in three unknowns 
has been reduced to two Poisson equations and one Helmholtz equation to 
be solved at each time step. Our first step concerning the solution 
procedure should involve solving the <j) equation (5-40) for a new value of 
<])*. The second step should be then to solve equation (5-41) for (<{)* + % t ) 
and after substitution for % t . At last we solve equation (5-42) for y t . Our 
history variables are ((), u, and v and they are updated after each time step 
as 



* 



<t> (t + At) = 2<)> - $ (t - At), 


(5-44) 


u (t + At) = 2 At [(% t ) x - (V t ) y ] + u (t - At), 


(5-45) 


v (t + At) = 2 At [(x t ) y + (\|/ t )J + v (t - At). 


(5-46) 



In other words, numerical integration of the three forecast equations 
involves 

a. solve first a Helmholtz equation for <{>, 

b. solve a Poisson problem for \|/, and finally 

c. solve a Poisson problem for %. 
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The space discretization consists of expanding the dependent variables 
in terms of basis functions defined on a variable mesh, and then 
orthogonalizing the error to the basis using the Galerkin procedure 
described in Chapter II. An appropriate approximating function for the 
rectangular subdivision should be a bilinear function (f). In this case the 
forecast set of equations in Galerkin form become 




4». f. 
y j j 

<D (At)" 



if i = J^ [(vQ) i f i4 <K j f j )] ) f i 



• J { £ [(uQ) j f j + £ (k j 9 ] ) + • Vj fj + £v» - v, 9 f. 



dy 



O At 



dy 






(5-47) 



O (At) 



/ v f * = - <uQ) j f 4 f i - J [ £r f j ] f i . 



(5-48) 



J v 2 [ £ (*)) f i + "’1 9 f i = J t £ [ (v® j f j - $r>j 9 ) f i 




(5-49) 
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The integral sign here means an area integral over the domain, the j 
subscript refers to Einstein summation for the dependent variables and the i 
subscript is the ith nodal equation. Following the integration-by-parts 
procedure the final form of the forecast equations is given by 






f. f. 

J J 



dy j dy 



O (At) 









-J{£ [(uQ) - f ; + K , ( f y ^ + <£>, 



dy 



O At 



[ v«t> - vij <^>j ) f i • J 777 <• • At > f i f i + J ir [u j (t ' A,) 



O (At) 



+ v j ( t -A, )< | ) j ] f i -jf 0 u j fj f i , 



( 5 - 50 ) 



J [<X>J <£>J + ( TT>> - - 1 t<"®j <§>; 1 f i 



' J [<vQ >j (57V f i • 



( 5 - 51 ) 



(5-52) 



S t Jr [<vQ) j f i • K i ( #y ) f i • 1 { lr [(uQ) j f i + K i } f i ■ 



The line integral along the north and south walls has been dropped 
from the vorticity equation (5-51), since the value of 9\y/9t is zero on the 
boundaries (\\f = constant). Also, the line integral along the north and south 
walls has been dropped from the divergence equation (5-52), because the 
value of the normal derivative along the north/south boundaries is zero. As 
we mentioned before the initial conditions should also satisfy the condition 
that the normal derivative of % along the north/south boundaries is zero. 



C. STABILITY ANALYSIS 

Here we analyze the primitive form of the forecast equations and a 
semi-implicit time scheme, since the results will be identical for the 
vorticity-divergence form, (Hinsman, 1983). The one-dimensional 
equations with a mean flow, U, are given by 



9u 9<{> 

aT + a^ 



= - u 



3u e 

5T +fv ’ 



(5-53) 



5v 

di 



dv 

3 ‘ u 3x ' f u ’ 



(5-54) 



3t dx 3x 



( 5 - 55 ) 
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Evaluating the time derivatives with a centered time differencing, and 
averaging the other terms on the left-hand side between time levels (t + At), 



and (t - At), equations (5-53), (5-54), and (5-55) become 

u(x,t + At) - u(x,t - At) 1 <)>(x + Ax,t + At) - <|)(x - Ax,t + At) 

2 At + 2 2 Ax 

<j)(x + Ax,t - At) - (})(x - Ax,t - At) n _ T T r u(x + Ax,t) - u(x - Ax,t) n 
+ 2A^ ^ J " ' U [ 2~Ax 1 



+ f v(x,t), 



v(x,t + At) - v(x,t - At) v(x + Ax,t) - v(x - Ax,t) 

= - U t — ] - f u(x,t). 



2 At 



2 Ax 



<{>(x,t + At) - 4>(x,t - At) 0 u(x + Ax,t + At) - u(x - Ax,t + At) 
2 At 2 2 Ax 



+ u(x + Ax,t - At) - u(x - Ax,t - At)^ _ ^ ^<j)(x + Ax,t) - <J)(x - Ax,t)^ 



2 Ax 



2 Ax 



Assuming now a function, F, given as 



F(x,t) = F' exp[i(kx + cot)]. 



and substituting into (5-56), (5-57), and (5-58), we can get 



( 5 - 56 ) 
( 5 - 57 ) 

( 5 - 58 ) 
( 5 - 59 ) 
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( 5 - 60 ) 



TT {”T ■ [f 2 + (k') 2 0 2 O] } = 0, 



At A 2 

At 



where 



s = sin(coAt) + GO u At, 



(6-61) 



c = cos(coAt), 



(k’) = s in(-^— )• 
Ax 



(5-62) 
( 5 - 63 ) 



The roots of equation (5-60) are given by 



s = 0, 



(5-64) 



s 2 = (fAt) 2 + c 2 (kAt) 2 <D. 



( 5 - 65 ) 



Requiring co to be real, the roots of (5-65) yield the following stability 
criterion 



At < 



1 






( 5 - 66 ) 
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J£L-I2LIXIAL_£GXI2IIHMS. 



A. TOPOGRAPHIC ROSSBY WAVE 

Let us first consider the horizontal momentum equations (3-8), and (3- 
9). If we cross differentiate (3-8) with respect to y, and (3-9) with respect 
to x, we obtain 



9 2 u 9u 9u 9 2 u dv 3u 9 2 u dv _ 9 2 h 

3y dt + dy dx + U dy dx + dy dy + V ^ : 2 dy ® dy dx’ 



( 6 - 1 ) 



3 2 v 9u dv d\ dv dv d\ 
dx dt + dx dx + U 0 x 2 dx dy V dx dy 



. du _ a 2 h 

dx ^ dx dy 



( 6 - 2 ) 



If h is eliminated from (6-1), and (6-2), then 



d dv 9u d dv 9u d dv 9u 

9t 9x 9y + U 9x 3x 9y + V 3y 9x 9y 



. ,9u 9u x dv du . ,9u dv x 

f( 3^ + a^ ) ' ( a7"37 )( aT + 



(6-3) 



Using now the following definition: 
^ dv du 

^ =0) * = a^"a7’ 



(6-4) 
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where £ is the vertical component of the vorticity, equation (6-3) yields 




(6-5) 



If we now use equation (3-15), equation (6-4) can be written in the 
following form 



where f is assumed to be constant. 

At this stage, in order to find the proper initial conditions for the 
particular case of the topographic Rossby wave, we can work out a 
simplified theory for Rossby waves (Phillips, 1965). 

Let us use cartesian coordinates and simplify the H variation, since H 
is a function of y, as 



dC _ C + f dH 
dt “ H dt 



( 6 - 6 ) 



or 




(6-7) 



H = D (1 - sy ), 



( 6 - 8 ) 
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where D is a constant, and 



s is the slope of H in the y-direction. 
Expanding the (£ + f) / H term as 



£ + f_ C + f _ (C+f)(l +sy) _ C + f + C sy + fsy 

H ~ D (1 - sy ) D D 



(6-9) 



and considering the term ( C, sy) to be very small, equation (6-7) becomes 



dt 



(C + f sy ) = 0. 



( 6 - 10 ) 



Equation (6-10) can be also written as 



dt dt. 



( 6 - 11 ) 



or 



^ + f sv = 0. 
dt 



( 6 - 12 ) 



Assuming small amplitude motion, and no mean flow, equation (6-12) 
becomes 
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+ f sv = 0. 
3t 



(6-13) 



In the case of a small Rossby number, we can also use the following 
relationships: 



C = V 2 y 



(6-14) 



and 



v 



8y 

37 ’ 



(6-15) 



where y = P / (pf). Equation (6-13) then becomes 



|<vSo + f,«§E)-0 



(6-16) 



or 




V) + P(^)=°, 



(6-17) 



where (3 = f s. 

Assuming a solution of the form 
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V = v(y) exp[ifi(x - ct)]. 



(6-18) 



substitution into equation (6-18) yields the following problem for \\f( y) 



,2 

,dw 2s 

( - (i ) (- 1 |!C) exp[ip(x - ct)] + P V (i|i) exp[i|i(x - ct)] = 0 

dy 


(6-19) 



+|-)\ir=0. 
dy 2 c 


( 6 - 20 ) 



Nondimensionalizing the domain of integration (as shown in Fig. 6.1), 
our boundary problem for y(y) becomes 



dV, 2 Pi 

dy 

with 


(6-21) 


< 

o 

II 

p 


( 6 - 22 ) 


V t (a) = V 2 (a), 


( 6 - 23 ) 
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Fig. 6.1 Schematic representation of the non-dimensional 
domain of integration. 
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(6-24) 



dVj(a) _ dy 2 (a) 
dy dy 



for the lower part of the channel, and 

d V, 2 fL 

— + =0, (6-25) 

dy 



with 



y 2 (i) = o, 


(6-26) 


y 2 (a) = ^(a). 


( 6 - 27 ) 


d\\f 1 (a) dy 2 (a) 
dy dy 


( 6 - 28 ) 


for the upper part of the channel. 


Here |ii, and P 2 are the correspondin 


wavenumbers. Also 




p! = s, f. 


( 6 - 29 ) 


P 2 = s 2 f. 


( 6 - 30 ) 
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where S\ is the slope of H in the y-direction for the lower part, and 
S 2 is the slope of H in the y-direction for the upper part. 



There are two distinct cases. In the first case, the lower part of the 
channel controls the phase speed c. The solutions to this case are 



Vj = A sin^y). 


(6-31) 


for the lower part, and 




Vj = B sinh[X 2 (1 - y)]. 


(6-32) 


for the upper part. For the particular case where a = 0.5 




Pi - - p 2 - 


(6-33) 


2 2 Pi 


(6-34) 


- (Xj) 2 = - 2n 2 - ^ . 


( 6 - 35 ) 



Furthermore, the coefficients A, and B are to be determined. In the second 
distinct case, the upper part controls the phase speed c. The solutions to 
this case are 
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V 2 = c sin[\ 2 (l - y)]. 



( 6 - 36 ) 



for the upper part, and 

Y 2 = D sinhC^y), ( 6 - 37 ) 

for the lower part. For the particular case where a = 0.5 




(6-38) 




(6-39) 




( 6 - 40 ) 



Once more, the coefficients C and D are to be determined. 

Let us now consider the second case, where the upper part controls the 
phase speed c. Boundary conditions (6-23), and (6-27) yield 

sin[X 2 (l - a)] C - [exp(X t a) - exp(- a)] D = 0. ( 6 - 41 ) 

In the same way, the boundary conditions (6-24) and (6-28) also yield 
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- X 2 cos[X 2 (l - a)] C - Xj [exp(Xj a) + exp(- Xj a)] D. 



( 6 - 42 ) 



Equations (6-41) and (6-42) lead us to the following eigenvalue problem 
Xj [exp(X t a) + exp(- X } a)] sin[X 2 (l - a)] 

+ X 2 [exp(Xj a) - exp(- X t a)] cos[X 2 (l - a)] = 0. ( 6 - 43 ) 

Using relationship (6-40), equation (6-43) becomes 

* * * *991 r> 

\ [exp(X 1 a) + exp(- Xj a)] sinf^Xj/ - ^ d * a)) 

* 9 9 1/9 * * 

+ [(Xj) - 2ii 2 ] [exp(Xj a) - exp(- Xj a)] 

* 9 9 1/9 

cos([(X i r - 2^] (1 - a)) = 0 , ( 6 - 44 ) 

where 

H 2 =~ W. (6-45) 

Here L is the horizontal length scale, and W represents the width of the 
channel. For our case we choose a channel 30° x 30° longitude by latitude, 
which gives 
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L = 5653510.75 m. 



W = 4896083.93 m, 
so |i2 equals to 5.4413981 

We are able to solve (6-44) numerically (see Table IV) using Newton's 
method. Using Table IV, we can determine coefficients C, and D for any 
particular mode. For instance if 

X* = 9.3191, 

X 2 = 5.2562, 

then in order to determine the relation between C, and D, the following set 
of equations has to be solved 

0.4912 C- 105.5791 D=0, (6-46) 

4.5787 C -984.0783 D =0, (6-47) 

resulting in 
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X o 


* 

X j 


5.2562 


9.3191 


11.1881 


13.5791 


17.3683 


18.9967 


23.6124 


24.8347 


29.8772 


30.8523 



Table IV Numerical solutions of equation (6-44), obtained using 
Newton’ s method. 
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( 6 - 47 ) 



0.4912 

105.5791 



C 



or 



C = 1.0, 



D = 0.004652 . 



Here C is chosen to be 1. Finally, the solution to the case where the upper 
part controls c, should be given by 

V 2 = sin [ 5.2562 (1 - y) ], 0.5<y<1.0, (6-48) 

for the upper part and 

\f 2 = 0.004652 [ exp ( 9.3191 y ) - exp ( - 9.3191 y ) ], 

0.0 < y < 0.5, ( 6 - 49 ) 

for the lower part. The graph of this solution is shown in Fig. 6.2. 

Following the same approach, with the use of Table IV, we are able to 
obtain 
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g. 6.2 Graph of the solution V 2 (y) where the upper part of 
the channel controls c, obtained from equations (6-48), 
and (6-49). 
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\y 2 (upper) = sin [ 1 1.1881 (1 -y) ], 0.5 < y < 1.0, 



( 6 - 50 ) 



Y 2 (lower) = - 0.000716 [ exp ( 13.5791 y ) - exp ( - 13.5791 y ) ], 

0.0 <y <0.5, 



for the case where 



\ = 13.5791, 
X 2 = 11.1881, 



shown in Fig. 6.3, and 



\|# 2 (upper) = sin f 17.3683 (1 - y) ], 0.5 < y < 1.0, 

(lower) = 0.0000506 [ exp ( 18.9967 y ) - exp ( - 18.9967 y ) ], 

0.0 <y <0.5, 



for the case where 



\ = 18.9967, 



(6-51) 



(6-52) 



(6-53) 
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Fig. 6.3 As in Fig. 6.2, except for equations (6-50), and (6- 
51). 
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x 2 = 17.3683, 
shown in Fig. 6. 4. 

The corresponding solutions for the case where the lower part controls 
the phase speed c, are given below 

Yj (lower) = sin ( 5.2562 y ), 0.0 < y < 0.5, ( 6 - 54 ) 

Yj (upper) = 0.004652 {exp[ 9.3191(1 - y)] - exp[ - 9.3191(1 - y)]}, 

0.5 < y < 1.0, (6-55) 

for the case where 
X\ = 9.3191, 

Xj = 5.2562, 
shown in Fig. 6. 5, 

Yj (lower) = sin ( 1 1.1881 y ), 0.0 <y <0.5, (6-56) 
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Fig. 6.4 As in Fig. 6.2, except for equations (6-52), and (6- 
53). 



103 






— _//=M. 

i*i i 




Fig. 6.5 Graph of the solution ^(y) where the lower part of 
the channel controls c. obtained from equations (6-54), 
and (6-55). 
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Vj (upper) = -0.0007 16{exp[l 3.579 1(1 - y)] - exp[ -13.5791(1 - y)]}, 

0.5 < y < 1.0, (6-57) 

for the case where 
X 2 = 13.5791, 

X 1 = 11.1881, 

shown in Fig. 6.6, and 

xj/^ (lower) = sin ( 17.3683 y ), 0.0 <y <0.5, (6-58) 

\\r l (upper) = 0.0000506{exp[18.9967(l - y)] - exp[ -18.9967(1 - y)]}, 

0.5 < y < 1.0, (6-59) 

for the case where 
X 2 = 18.9967, 

X x = 17.3683, 
shown in Fig. 6. 7. 
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Fig. 6.6 As in Fig. 6.5, except for equations (6-56), and (6- 
57). 
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Fig. 6.7 As in Fig. 6.5, except for equations (6-58), and (6- 
59). 
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In our experiments we select initial conditions in the form described 
by equations (6-54) and (6-55). That choice appears to be an excellent one, 
providing us with the most proper modes to observe the topographic 
Rossby wave. It is desired that the initial conditions allow a relative 
amount of control for the input parameters as well as satisfing the boundary 
conditions. 

The next logical step is to determine the analytic expression for the 
streamfunction \| /. Following closely equations (6-54) and (6-55), the exact 
expression for \| /, is given by 



v = 4" t sin (— y)] t sin (- 

2 w 



2k n 



L 




( 6 - 60 ) 



for the lower part (0.0 < y < 0.5), and 



A. * * 2tc n 

\\f = 0.004652 — [exp(X 2 y) - exp( - X 2 y)] [sin(— x)] 




(6-61) 



for the upper part (0.5 < y < 1.0), where 



A = amplitude of perturbation, 

W = width of the channel (4896083.93m), 
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L = length of the channel (5653510.75m), 
n = wave number, 

U m = mean flow speed, 

y m id - middle point of the channel, 

?ii = 5.2562, and 
X 2 * = 9.3191 

The first term in the expressions (6-60) and (6-61) represents the 
perturbation part. The second term is the north/south slope necessary to 
support a mean flow of U m - The last term plays the role of the mean depth 
term. The geopotential height, 0, is related geostrophically to the 
streamfunction, \j/, by the following relationship 




( 6 - 62 ) 



resulting in 




( 6 - 63 ) 



for the lower part (0.0 < y < 0.5), and 
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4> = 0.004652 [expCa^y) - exp(- a^y)] sin(a 2 x) 



- f 0 U m (y 



y^ + <t> ' 



(6-64) 



for the upper part (0.5 < y < 1.0), where 



a 



l 




( 6 - 65 ) 



a 



2 



2k n 



( 6 - 66 ) 




(6-67) 



The u, and v components of velocity can be derived using the 
following geostrophic expressions 



u = - 



1 3<J> 

7 W’ 



( 6 - 68 ) 



v = 



1 3<j) 
f 9x 



( 6 - 69 ) 



resulting in 
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(6-70) 



A a 

u (lower) = U m — - — cos^y) sin(a 2 x), 0.0 < y < 0.5, 
A a 2 

v (lower) = — ^ sin^y) cos(a 2 x), 0.0 < y < 0.5, 



and 



u (upper) = U m - 0.043352 — [exp(a 3 y) + exp(- a 3 y)] sin(a 2 x), 

0.5 <y < 1.0, 



A a- 

v (upper) = —— [exp(a 3 y) - exp(- a 3 y)] cos(a 2 x), 0.5 < y < 1.0. 



The initial vorticity can be also derived, using the 
relationship 



C-vV 



resulting in 




[(a t ) 2 + (a^ 2 ] sin(ajy) sin(a 2 x), 



for the lower part (0.0 < y < 0.5), and 



(6-71) 

(6-72) 

(6-73) 

following 

(6-74) 

(6-75) 
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C = 0.004652 y [(a^ 2 - (c^) 2 ] [exp^y) - exp(- ^y)] sin^x). 



(6-76) 



for the upper part (0.5 < y < 1.0). 

We set the initial divergence equal to zero, assuming the initial fields 
to be almost geostrophic. The initial fields of geopotential, u-component, 
v-component, vorticity, and divergence for the case where the lower part of 
the channel controls the phase speed, c, and no topography effect is 
involved (Equations 6-60 through 6-76), are illustrated in Fig. 6.8. 

B. HYDRAULIC JUMPS 

The stationary theory predicts the formation of jumps in pairs, one 
upstream and one downstream of the rigid ridge (Houghton and Kasahara, 
1968). The classification shown in Fig. 4.5, is strictly valid for non- 
rotating flows. No equivalent theory exists for flows over mountains in a 
rotating system. Houghton (1969) and Williams and Hori (1970) consider 
the transient motion of a shallow water layer on an f-plane without 
mountains starting from an initial velocity disturbance of magnitude U over 
a length L. 

For our case we consider a nonrotating Shallow-Water system (f = 0). 
It is desired once more, that the initial conditions allow a relative amount 
of control for the input parameters as well as satisfing the boundary 
conditions. The forecast model history-carrying variables are <J), u, and v. 
The analytic expression for the sreamfunction, y, is given by 
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Fig. 6.8 Initial conditions for the GFEM model with a 
rectangular subdivision, and wave number one. 
Contour intervals are 600 m 2 /s 2 for geopotential 

height, 0.2 m/s for u and v, 0.6 x 10‘ 6 S ’ 1 for 
vorticity. Nodal points are denoted by an x. 
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y = - U y, 0.0 < y < W, 



( 6 - 77 ) 



where U is the velocity of the flow in the region far from the ridge. Also 
the initial geopotential height, 4>o, is set equal to the mean depth, gH, 

Finally, the initial u- and v- components of the velocity are given by 



= u, 


(6-78) 


= 0. 


(6-79) 
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XII t .EXEEKIMEKIS.-A£LJ2.KES.ULIS. 



Our first experiment involves bottom topography which is composed 
of two regions. These two regions can provide us with either an east-west 
oriented ridge or valley. We consider conditions with no mean flow, the 
objective being to examine how well our model simulates the topographic 
Rossby wave, by comparison with the theoretical phase speed values. 

Our second experiment is to investigate the ability of the same model 
to create hydraulic jumps analogous to that predicted from the analytical 
approach (Chapter IV). In this case we consider a mean flow forced to 
pass over a topographic ridge which extends north-south across the 
channel. In this experiment we will consider several distinct cases 
corresponding to different discrete domains obtained by the theory (Fig. 
4 . 5 ). 



A. EXPERIMENT I 

We perform experiment I using the GFEM rectangular model described 
in Chapter VI. The basic difference between the rectangular and triangular 
models is in the approximating polynomials. The rectangular polynomials 
are bilinear while the triangular polynomials are linear. Many integrals 
require evaluation during the integration process. We could use numerical 
quadrature, however, a more efficient method is available through the use 
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of natural coordinates. We can accomplish quadrature with no error by 
formula with this method. A description of the natural coordinate method is 
given in Appendix A for the rectangular discretization. 

The diagram for the grids for the rectangular subdivision is shown in 
Fig. 5.1. The domain of integration is 5,653.5 km in the x-direction and 
4,896.1 km in the y-direction. The model has 12 increments in the x- and 

k 44 

y- directions, which gives the model 156 degrees of freedom. The Ax is 
471. 1 km, and the Ay is 408.0 km. The value of the Coriolis parameter is 
taken to be 0.00010284, corresponding to 45.0* N latitude. 

The initial conditions for the case of the topographic Rossby wave are 
described in Chapter VI. A small wave perturbation is added to the 
geopotential field which includes the mean height and the required 
north/south slope in case of non zero mean flow. It consists of a wave with 
a wavenumber one, confined primarily to the south part of the channel 
domain. In our case, we choose a mean depth of 1,000 meters, and no 
mean flow. The motion is confined in a channel with cyclic boundary 
conditions as shown in Fig. 7.1. We examine small amplitude wave motion 
encountering different slopes of the bottom topography in the y-direction as 
illustrated in Fig. 7.2, and 7.3. We can visualize the whole setting as if we 
placed a long triangular mountain with its peak centered in the middle of 
the width of the channel. It is obvious that the slope of the lower part of 
the channel, s\, corresponds to the south slope of the mountain, while the 
slope of the upper part, S 2 , corresponds to the north slope. No matter what 
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Fig. 7. 1 Schematic representation of the domain of integration. 
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n "i 




Fig. 7.2 Schematic representation of the bottom topography for 
the case of triangular mountain. 
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Fig. 7.3 As in Fig. 7.2, except for the case of triangular valley. 
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the case is, the following relationship holds: 



Sj = - s 2 , ( 7 - 1 ) 

where Si is assumed to be greater than or equal to zero for most of 
the cases. 

B. RESULTS I 

We integrate the forecast equations over a time interval of 96 hours, 
and we plot the results every 48 hours. Our first concern is to examine the 
case of no topography. In this particular case I, the peak of the triangular 
mountain is zero, and both slopes, si and S 2 , are equal to zero as we can 
see from Table V. The initial conditions are almost purely geostrophic as 
clearly illustrated in Fig. 7.4. The integration produces forecast fields 
almost identical to the initial fields, as shown in Figs. 7.5, and 7.6. That is 
expected, since no topographic or beta effect is involved. If we now 
increase slightly the magnitude both of S\ and S 2 , as shown in Table V for 
case II, the 48 hour integration does not show any significant change in the 
forecast fields, as we can see in Fig. 7.7. However, the 96 hour integration 
yields a very small tendency for a westward shifting valid for all the 
forecast fields, and for the lower part of the channel. 

Our next step, case III, is to increase the peak of the mountain to 
163. 1 m. In both the 48, and 96 hour integrations, a tendency for the lower 
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Peak (6) 


Peak (z) 


Slope (south) 


case I 


0.0 


0.0 


0.0 


case II 


400.0 


40.8 


0.0000167 


case IE 


1600.0 


163.1 


0.0000666 


case IV 


4800.0 


489.3 


0.0001999 


case V 


- 400.0 


-40.8 


-0.0000167 


case VI 


- 1600.0 


- 163.1 


- 0.0000666 


case VC 


- 4800.0 


- 489.3 


-0.0001999 



Table V Peak (in geopotential meters and in meters), and south 
slope values for cases I through VII. 
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Fig. 7.4 Initial conditions for experiment I (cases I through 
VII). Contour intervals are 600 m 2 /s 2 for geopotential 
height, and 0.2 m/s for u and v. Nodal points are 
denoted by an x. 
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Fig. 7.5 



As in Fig. 7.4, except for case I, 
integration. 



and a 48 hour 
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Fig. 7.6 As in Fig. 7.5, except for a 96 hour integration. 
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Fig. 7.7 



As in Fig. 7.5, except for case 



II. 
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Fig. 7.8 As in Fig. 7.6, except for case II. 
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part forecast fields for a westward shifting is clearly noticed, as well as a 
small tendency for eastward shifting corresponding to the upper part fields 
(Figs. 7.9, and 7.10). In case IV, which is our last case of positive south 
slope si, we increase the peak of the mountain to 489.3 m. This value 
corresponds to almost the half of the mean depth value considered for our 
shallow water approximations. In both the 48, and 96 hour integration it is 
very evident the significant westward sifting of the forecast fields for the 
the lower part of the channel, and the eastward shifting of the same fields, 
characterizing the upper part, as we can see in Figs. 7.11, and 7.12. The 
results obtained here, will be compared quantitatively with the 
corresponding analytical values. 

At this stage we wish to examine cases V, VI, and VII, all having 
negative south slopes. The above mentioned cases correspond to cases II, 
III, and IV, but with exactly opposite sign slopes. We can visualize the 
whole setting here as if we placed reverse triangular mountains (valleys) of 
different heights with their lowest points centered in the middle of the 
width of the channel, and their peak values always to be given by negative 
values. The 48 hour integration results for case V do not show any 
significant change in the forecast fields, as we can see in Fig. 7. 13, while 
the 96 hour integration yield a very small tendency for a eastward shifting 
valid for all the forecast fields, and for the lower part of the channel (Fig. 
7.14). In case VI, the valley is 163.1 m deep in the center. In both the 
48, and 96 hour integrations, a tendency in the lower part forecast fields 
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As in Fig. 7.5, except for case III. 
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Fig. 7. 10 As in Fig. 7.6, except for case III. 
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Fig. 7.11 



As in Fig. 7.5, except for case IV. 
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Fig. 7. 12 As in Fig. 7.6, except for case IV. 
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Fig, 7. 13 As in Fig. 7.5, except for case V. 
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As in Fig. 7.6. except for case V. 
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for a eastward shifting is evident, as well as a small tendency for westward 
shifting corresponding to the upper part fields (Figs. 7.15, and 7.16). The 
same tendency becomes more evident in both the 48, and 96 hour 
integration, i. e., the significant westward sifting of the forecast fields for 
the the lower part of the channel, and the eastward shifting of the 
corresponding upper fields, as it is clearly shown in Figs. 7. 17, and 7. 18. 

At this point, we strongly believe that there is a certain link between 
the presence of topography and the observed shifting in all the forecast 
fields, because no shifting at all is been observed in the flat case of no 
topography. That specific link has to be the topographic Rossby wave 
whose existence requires the topographic effect in a rotating system (with 
constant f ). The reason is that the Rossby wave in general, can exist only 
in the presence of an ambient potential-vorticity gradient. In case I, of 
course, no potential-vorticity gradient is present, and that is why we do not 
observe any sign of the Rossby wave. In case II, III, and IV, the positive 
y-direction (shown in Fig. 7.19), is also the direction of increasing 
ambient potential-vorticity. If we consider a fluid column initially at rest, 
but later to be displaced in the positive y-direction, then in order to 
conserve its total potential-vorticity, the wave potential-vorticity C,o ' Fr[o> 
has to be decreased. This will balance the excess of the ambient potential 
vorticity in its new place. There are two ways for this to be accomplished. 

a. The fluid will have a tendency to be squeezed, since its new 
position appears to be shallower than the previous one. This will induce the 
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Fig. 7. 15 As in Fig. 7.5, except for case VI. 




Fig. 7. 16 As in Fig. 7.6, except for case VI. 
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As in Fig. 7.5, except for case VII. 
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Fig. 7. 18 As in Fig. 7.6, except for case VII. 
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Fig. 7.19 The required ambient potential-vorticity gradient. A 
clue for the physical explanation of the topographic 
Rossby wave oscillation. 



139 



production of negative vorticity due to the vortex-tube compression. 

b. The act of the squeezing can not be complete, since the upper 
surface is not restricted or bounded. In this case the column will have a 
tendency to ride at least partially up, the slope resulting in a greater value 
of TJo in its new position than its neighbors. 

It is important here to note that both effects, Co < 0. and rjo > 0, act to 
reduce the quantity Co * Fr|o- Both effects also, will result in a clockwise 
circulation in the fluid around the column. The clockwise circulation in the 
fluid column C, will force the adjacent column to its right, R, into deeper 
fluid, and the adjacent column to the left, L, to be squeezed into shallower 
fluid, as shown in Fig. 7. 20. The column R, will become the center of a 
counter-clockwise circulation, while column L, will become the center of a 
clockwise circulation. Both contributions from columns L, and R, will 
force the return of column C, toward its original position. This will result 
in an overshoot due to the column C inertia, and the oscillation will 
continue. This is a very simplified view of the phenomenon but it clearly 
shows a very important aspect. The strength of the restoring mechanism 
depends on the vigor of the circulation induced on neighboring fluid 
columns by the displaced column. Following the same approach, in cases 
V, VI, and VII, the positive y-direction is the direction of decreasing 
ambient potential-vorticity, and everything described above applies in the 
opposite sense. 



140 




Fig. 7.20 The position of the three-point vortices L, C. and R 
at three successive times. Initially collinear and 
positioned along an isobath, C is displaced 
upwards, producing velocities at L and R which 
move them as shown. The vorticity induced on L 
and R produces a velocity at C with a tendency to 
restore it to its original position. 
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From a purely theoretical point of view, if we recall equation (3-47) we can 
easily find out that for positive values of slope s the phase speed in the x- 
direction is always negative or, referring to our case, westwards. For the 
opposite case of negative values of s, the propagation of the Rossby wave 
in the x-direction is always positive or eastwards. Also, from the same 
equation (3-47), it is obvious that for increasing values of s (in 
magnitude), the propagation phase speed also increases. In other words, if 
we keep increasing the peak of our triangular mountain we should expect 
the presence of the topographic Rossby wave to become more and more 
evident. In order to examine the actual phase speed in more detail, we can 
Fourier analyze the v-component field, and obtain the wave component 
phase speed every 24 hours. The observed phase speed is then given by 




L(y<t> 2 ) 

2 k (tj - 1 2 ) 



(7-2) 



where L is the length of the channel. These phase speeds, averaged over 96 
hours for cases II through VII, are given in Table VI, and in Figs. 7.21, 
and 7.22. 

The analytic phase speed given by (3.41) can be rewritten in the 
following form: 
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Peak (z) 


C (estimated) 


C (corr. fs) 


C (corr. fs+H) 


C (observed) 


case II 


40.8 


-0.72 


-0.49 


-0.49 


-0.49 


case m 


163.1 


-2.87 


- 1.97 


- 1.98 


-1.99 


case IV 


489.3 


-8.61 


-6.14 


-7.12 


-7.44 


case V 


-40.8 


0.72 


0.49 


0.49 


0.49 


case VI 


- 163.1 


2.87 


1.97 


1.85 


1.74 


case VTI 


-489.3 


8.61 


6.14 


5.22 


4.89 



Table VI Estimated, estimated corrected by the free surface term, 
estimated corrected by both the free surface and H 
terms, and observed phase speed values (in m/s), for 
cases II through VII. 
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Fig. 7.21 Comparison of the observed phase speed values, 
with estimated, estimated corrected by the free 
surface term, and estimated corrected by both the 
free surface and H terms, phase speed values (in 
m/s) for cases II through IV (scatter diagram). 
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As in Fig. 7.21, except for cases V through VII. 
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(7-3) 



where hM is the mountain height, and the term [(|ii) 2 + (Xi) 2 ] is non- 
dimensional. The analytic phase speeds resulting from (7-3) are also given 
in Table VI. 

In general, the analytic phase speeds are all larger in magnitude than 
the corresponding model phase speeds, obtained from (7-2). The reason 
why this happens is that the analytic theory developed here is based on a 
rigid upper lid, but at the same time, the GFEM model used for our 
forecasts has an upper free surface. If we wish to include the upper free 
surface effect in equation (7-3), the term [(fo) 2 / gH] must be added to the 
denominator, resulting in 



f l h M 

° H W 



C = 



2 



‘ (7-4) 




w 2 gH 



where H represents the mean depth value. 
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The value of the new added term is about the half of the denominator 
value, so that it will reduce the analytic phase speed value by around 30%, 
which would bring C and Cp, into more agreement for most of the cases. 
However, the highest mountain peak case, case IV, or the lowest depth 
valley case, case VII, can not be satisfied by only considering this change. 
The above mentioned agreement could be improved for case IV, by using a 
smaller H value into (7-4), since the true average depth in this case is less 
than 1 km. Using the same arguments for case VII, the improvement of the 
desired agreement could be accomplished by using a larger H value into (7- 
4), since the true average depth in this last case is greater than 1 km. 

Table VI and Figs. 7.21, 7.22, 7.23, and 7.24, compare the observed 
phase speeds with the various theoretical estimates. As expected, each new 
correction improves the agreement with the model phase speed. The phase 
speed for the shallower topography are extremely accurate. The reason for 
the lack of agreement for larger H values is due to the uncertainty for what 
the correct value of H is to use in equation (7-4). Therefore, for higher 
values of the bottom topography the error is larger because using the 
correct value of H in equation (7-4) is most important. Clearly, experiment 
I shows that our numerical model is able to handle the topographic Rossby 
wave case extremely well. 
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As in Fig. 



7.21, except for bar representation. 
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Fig. 7.24 As in Fig. 7. 22, except for bar representation. 
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C. EXPERIMENT II 

We perform experiment II using the same GFEM rectangular model, as 
in experiment I, the only difference is in the resolution of the model. The 
model now has 24 increments in the x- and y- directions which gives the 
model 576 degrees of freedom. The domain of integration is 727.6 km in 
the x-direction, and 630. 1 km in the y-direction. The value of Coriolis 
parameter is taken to be zero for all of the cases, corresponding to a 
nonrotating system. 

The initial conditions for experiment II are described in Chapter VI. 
The boundary conditions in x are chosen to be periodic once more, that is 



These boundaries at x = 0, and x = L, are placed sufficiently far from the 
ridge so that the desired asymptotic conditions are well established in the 
vicinity of the ridge before wave motions are able to be fed back into this 
region by the periodic boundary conditions. The height profile of the 
orographic ridge is given by 



u (0) = u (L), 
<t> (0)=<}> (L). 



(7-5) 

(7-6) 




(7-7) 
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where Z is the width and Hm the height of the mountain, as shown in Fig. 
7.25. 

D. RESULTS II 

We integrate the forecast equations over a maximum time interval of 
50.6 minutes and we plot the results at t = 5.6, 16.9, 28. 1, 39.4, and 50.6 
minutes. The five distinct cases we run (I through V), are summarized in 
Tables VII and VIII. Cases III and IV are expected to produce a hydraulic 
jump because they lie in domain II (see Fig. 4.5). In each case the 
equations are integrated from an initial state where u and h are both 
uniform. The u-field for case I is shown in Fig. 7.26. It clearly shows the 
rapid development of a speed maximum over the center area of the mountain 
and slightly on the lee side. A secondary speed maximum also forms over 
the ridge area and it moves upstream with time. The 0- field , shown in Fig. 

7.27, indicates the earlier development of low pressure-field over the lee 
side of the ridge. The high pressure-field over the east part of this pattern 
has an obvious tendency to move upstream with time. Since the main 
feature of the flow, the speed maximum over the center area of the 
mountain, is quite stable we regard this as a no-jump case in agreement 
with the theory presented before in chapter IV. The u-field for case II, a 
case with slightly higher mountain, and stronger mean flow, shown in Fig. 

7.28, is generally similar to case I except that the perturbation now 
represents a larger fraction of the mean flow. The same arguments appear 
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Fig. 7.25 Schematic representation of the position of bottom 
topography along x-axis, valid for each node per 
horizontal row, for the hydraulic jump case 
(experiment II). 
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Fo 


R 


F 


H (m) 


U (m/s) 


D 


case I 


0.2 


0.1 


1.010 


1000.0 


19.80 


I 


case 13 


0.4 


0.2 


1.040 


1000.0 


39.60 


I 


case HI 


0.3 


0.7 


1.220 


1000.0 


29.69 


n 


case IV 


0.8 


0.2 


1.149 


500.0 


56.00 


n 


case V 


1.4 


0.05 


1.407 


400.0 


87.65 


in 



Table VII Froude number (Fo), maximum height of the ridge (R), 
parameter F, mean depth (H), mean flow (U), and 
domain (D), for cases 1 through V. 
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F 


R 


D 


Jump case ? 


case I 


1.010 


0.10 


I 


no 


case II 


1.040 


0.20 


I 


no 


case m 


1.220 


0.70 


II 


yes 


case IV 


1.149 


0.20 


II 


ves 


case V 


1.407 


0.05 


III 


no 



Table VIII Parameter F, maximum height of the ridge (R), 
domain (D), and classification of the asymptotic 
flow conditions, for cases I through V. 
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.26 U-component amplitude as a function of time for case I. 
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.27 O-amplitude as a function of time for case 1. 
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Fig. 7.28 As in Fig. 7.26, except for case 11. 



to be valid for the 0-field also, shown in Fig. 7.29. We also regard case II 
as a no-jump case in agreement with the theory. 

The u-field for case III is shown in Fig. 7.30. In this case the wind 
maximum on the lee side of the mountain, continues to grow with the time. 
At the same time the 0-field, shown in Fig. 7.31, indicates that the low 
pressure-field centered on the lee side, changes rapidly to nearly zero. We 
regard case III as a jump case because it is not approaching steady state. 
The resolution of the model is too poor to allow a detailed description of 
the small- scale jump zone. The theory also indicates that this is a jump 
case. Fig. 7.32 contains the u-field for case IV. The 0-field for the same 
case is given by Fig. 7.33. The behavior of case IV is similar to case III, 
so that, this is also a jump case in agreement with the theory. 

The fields of u and 0 for case V, are shown in Figs. 7.34 and 7.35 
respectively. These field patterns are similar to those of cases I and II, the 
only difference is a downstream shifting for both the u- and 0- fields. The 
growth in the amplitude of the curves appears to be stabilized, so that we 
regard case V as a no-jump case in agreement once more with the theory. 

In all the investigated cases (I through V), jump formation is indicated 
by both u and 0 amplitudes which continue to increase with time. In order 
to study the behavior of each one of the jump cases in more detail, better 
resolution will be required as well as a larger domain to reduce the 
boundary effects. 
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Fig. 7.34 As in Fig. 7.26, excepi for case V. 
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XII£.££U.£LIIS.I£IiS. 



In this thesis we have investigated how well a particular shallow-water 
finite element prediction model handles surface topography. In both 
experiments the flow is confined within an east-west channel with periodic 
boundary conditions. 

In the first experiment the bottom is composed of two regions of 
constant and opposite north-south slope with no east-west variation, so that 
the bottom is either an east-west ridge or an east-west valley. In order to 
have the proper initial conditions, the analytic Rossby wave solutions are 
derived. These are obtained by solving the linearized quasi-geostrophic 
equations with the free surface assumed to be rigid. Each solution has a 
sinusoidal variation with y over one bottom slope and exponential decay 
over the other slope. With the simple Rossby formula the direction of 
propagation is determined by the bottom slope in the region which has the 
largest wave amplitude. These solutions are similar to the trench wave 
solutions obtained by Mysak et al (1979), who used piecewise exponential 
bottom profiles. The numerical solutions with the initial conditions given 
by the linear solutions produced smoothly propagating solutions. The phase 
speeds were very accurate when the Rossby formula was corrected for free- 
surface effects and mean depth. 

In the second experiment a north-south ridge was placed across the 
channel with sine-squared east-west variation. The Coriolis parameter was 
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set to zero and the initial u and <|> were constant. The theory of Houghton 
and Kasahara (1968) for the formation of hydraulic jumps was reviewed. 
The equations were integrated for five initial conditions. In each case a 
speed maximum occurred over or just downstream from the ridge and low 
heights were found on the lee side of the ridge. For three of the cases the 
solutions reached an approximate steady state. These agreed with the theory 
which predicted no jumps. The two other cases lead to increasing winds 
and decreasing heights with no steady state. These were jump cases 
according to the theory. In these two last cases the model resolution was 
inadequate to simulate the formation of the hydraulic jumps in detail. 

Finally, the finite element model performed well for two very different 
topographic effects. Further testing is required on the jump cases with 
higher resolution. Furthermore, the effect of the semi-implicit 
discretization on the hydraulic jumps should be determined. 
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APPENDIX 



NUMERICAL QUADRATURE 



A very fast and efficient method is required for the evaluation of the 
integrals obtained by the Galerkin approach. For the rectangular 
subdivision, integration formulas are based on an orthogonal axis 
transformation. In this case, the integrals to be evaluated contain either 
products of the basis functions, products of derivatives of basis functions, 
or a mixture of both. 

An orthogonal axis transformation will allow quadrature formulas for 
rectangles. We are able to transform the rectangle shown in Fig. A. 1 by 
using the following 







( A - 1 ) 



r| = 



y-y 0 

b * 



( A - 2 ) 



where the values of and T) at each corner are shown in parentheses. 
Using fj as a basis function, we can express f; as 

( 1 + C) ( 1 + Tl; Tl ) 

f.= . (A- 3) 

1 4 
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Fig. A. 1 Orthogonal axis transformation for rectangular 

integration formulas. 
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Derivatives of the basis functon are given by 



df. 1 df. 

i _ 1 1 

3x a 3 ^ 



and 



3f. 1 3f. 

1 __ 1 1 

3y “ ¥ ar, ' 



( A - 4 ) 



( A - 5 ) 



The interaction coefficients can now be determined for a derivative or 
straight inner product. The straight inner product is given by 

l l 

c^JJ fjfjdxdy-abJ j fjfjdCdn 

-1 -1 



1 1 

= ^ J (1 + C i C )(1 + C J C )d cJ ( 1 + T i i T i )( 1 +Ti j Ti )d n 

-1 -1 



76 (2 + j^i^j )(2 + y'ni'Hj >• 



( A - 6 ) 



Also, the mixed derivative is given by 
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9f. 9f. her or. at. 

^ dx dy = — — — dC dr\ 

ox dx a J J y 7*r 



1 1 



9f. 9f. 



-l -l 



ac ac 



i i 

= ^ J C; Cj d C J ( 1 + Tli ^1 ) ( 1 + Tlj T1 ) dT ! 

-1 -1 



= — r— (2^.^.)(2 + ~r\.r\.). 
16a 3 1 J 



( A - 7 ) 
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